home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / VideoToolbox 96.06.15 / VideoToolboxSources / ConvolveX.c < prev    next >
Text File  |  1996-01-25  |  12KB  |  294 lines
  1. /* 
  2. ConvolveX/Y.c
  3. © 1989-1993 Denis G. Pelli
  4. Simple one-dimensional convolution. Convolve image in srcBits by array f[] to
  5. produce image dstBits. src and dst may have any pixelSize, but the calculation
  6. itself is done with 8 bits per pixel. Assumes that everything has its origin at
  7. its center. Before calling ConvolveX/Y you should call SetEntries() to set up
  8. the color tables appropriately (unless source and destination pixmaps are both 8
  9. bits/pixel, in which case the color tables are ignored).
  10.  
  11. ConvolveX/Y fills each pixel in the destination with a weighted sum over the
  12. source, centered on the corresponding pixel in the source. To avoid edge
  13. effects, it only computes pixels in the destination which have full support in
  14. the source. If necessary, ConvolveX/Y will inset the source or destination rect,
  15. so that the destination width/height plus the array length equals the source
  16. width/height.
  17.  
  18. It would be nice to convert this to use Get and SetPixelsQuickly instead of
  19. directly accessing the pixels. This would introduce a small overhead, since the
  20. numbers would be copied into a long array and from thence into the pixels, instead
  21. of being stuffed directly, but would probably make this program quite a bit easier
  22. to read.
  23.  
  24. BUGS:
  25. Setting IMPROVED_ROUNDING to true improves results in some applications, but
  26. screws up others (e.g. the demo Filter.c). I don't have time to figure this out at 
  27. the moment, so I'm leaving a flag that can be turned on or off to suit your 
  28. application. I expect to fix this and eliminate the flag in the future.
  29.  
  30. I think it overflows if you try to filter 16 or 32-bit deep images. 
  31.  
  32. HISTORY:
  33. 4/1/89 dgp    twice as fast, thanks to recoding to use long instead of float arithmetic.
  34. 4/2/89 dgp    made universal by using line buffers for both src and dst, and interfacing
  35.             to the src and dst PixMaps solely through CopyBits.
  36. 4/2/89 dgp    ConvolveX/Y.c, a pre-processor conditional was added to generate ConvolveX
  37.             and ConvolveY from the SAME source. The two files are IDENTICAL except
  38.             for the preprocessor definition of CONVOLVE_X as 1 to generate ConvolveX,
  39.             or 0 to generate ConvolveY.
  40. 10/9/89 dgp    Fixed overflow bug. Speeded up the convolution loop slightly.
  41.             Added conditional to use CopyBitsQuickly instead of CopyBits when source and
  42.             destination pixelSize are 8 bits. This has two virtues. It's faster, and
  43.             it bypasses the Color Manager's transformations through the color and inverse
  44.             color tables, which costs time and loses accuracy.
  45. 10/10/89 dgp Now use RectToAddress to determine the pixelSize of src and dst, to deal
  46.             correctly with window's pixmap.
  47. 7/25/91    dgp    Added explicit preprocessor symbol to turn diagnostic printout on and off.
  48. 8/24/91    dgp    Made compatible with THINK C 5.0.
  49. 8/27/92    dgp    Added Gestalt().
  50. 11/10/92 dgp Now initialize all fields of pixmap, especially the PMVersion.
  51.             Identify all error messages as originating from here: ConvolveX/Y.
  52.             Replaced compile-time conditional by a run-time conditional, so
  53.             that the ConvolveY.c file is no longer needed.
  54. 1/22/93    dgp    SwapMMUMode().
  55. 2/1/93    dgp Improved rounding. Removed calls to Gestalt, as unnecessary overhead.
  56. 2/7/93    dgp    Introduced switch IMPROVED_ROUNDING to disable rounding for debugging.
  57. 2/8/93    dgp Fixed overflow problem caused by the new rounding.
  58. 7/9/93    dgp check for 32-bit addressing capability.
  59. 6/18/94    dgp    can32 is now computed by calling TrapAvailable(_SwapMMUMode), which 
  60.             returns the correct answer even on Macs with dirty ROMs.
  61. 5/23/95 dgp Apple changed the prototype in the header file from SwapMMUMode(char *) to 
  62.             SwapMMUMode(signed char *). To retain compatibility with both old and new
  63.             headers, I cast the argument (void *).
  64. */
  65. #include "VideoToolbox.h"
  66. #ifndef __TRAPS__
  67.     #include <Traps.h>        // _SwapMMUMode
  68. #endif
  69. void ConvolveXY(double f[],int dim,BitMap *srcBits
  70.     ,BitMap *dstBits,Rect *srcRectPtr,Rect *dstRectPtr,int convolveX);
  71. #define DIAGNOSTICS 0        // true or false
  72. #define IMPROVED_ROUNDING 1    // true or false
  73.  
  74. void ConvolveX(double f[],int dim,BitMap *srcBits
  75.     ,BitMap *dstBits,Rect *srcRectPtr,Rect *dstRectPtr)
  76. {
  77.     ConvolveXY(f,dim,srcBits,dstBits,srcRectPtr,dstRectPtr,1);
  78. }
  79.  
  80. void ConvolveY(double f[],int dim,BitMap *srcBits
  81.     ,BitMap *dstBits,Rect *srcRectPtr,Rect *dstRectPtr)
  82. {
  83.     ConvolveXY(f,dim,srcBits,dstBits,srcRectPtr,dstRectPtr,0);
  84. }
  85.  
  86. void ConvolveXY(double f[],int dim,BitMap *srcBits
  87.     ,BitMap *dstBits,Rect *srcRectPtr,Rect *dstRectPtr,int convolveX)
  88. {
  89.     Rect mysrcRect,mydstRect;
  90.     int dstLines,srcLines,dstLength,srcLength;
  91.     int lines;
  92.     unsigned char *dstPtr=NULL;
  93.     register unsigned char *srcPtr=NULL;
  94.     register int i;
  95.     Rect dstRect, srcRect;
  96.     register long t;
  97.     register long *ifun, *fPtr;
  98.     long ifunScale;
  99.     int x,y,srcx,imin,imax;
  100.     double tmax;
  101.     PixMap srcTmp,dstTmp;    /* line buffers */
  102.     #if DIAGNOSTICS
  103.         GDHandle device=NULL;
  104.         double a;
  105.     #endif
  106.     short srcPixelSize,dstPixelSize;
  107.     int error;
  108.     long value;
  109.     signed char mmuMode=true32b;
  110.     Boolean can32;
  111.     
  112.     value=0;
  113.     error=Gestalt(gestaltAddressingModeAttr,&value);
  114.     can32=TrapAvailable(_SwapMMUMode);
  115.     RectToAddress((PixMap *)srcBits,srcRectPtr,NULL,&srcPixelSize,NULL);
  116.     RectToAddress((PixMap *)dstBits,dstRectPtr,NULL,&dstPixelSize,NULL);
  117.     ifun = (long *) malloc(dim*sizeof(*ifun));
  118.     if(ifun == NULL) {
  119.         PrintfExit("ConvolveX/Y: Sorry, couldn't allocate enough memory for working array\007\n");
  120.     }
  121.     /*
  122.     We will sum in a long register.
  123.     The sum over the point spread function f[] will yield at most tmax, defined below,
  124.     so we can safely multiply f[] by LONG_MAX/tmax, without danger of overflowing t.
  125.     */
  126.     for (tmax=0.0, i=0;i<dim;i++) tmax += fabs(f[i]);
  127.     tmax *= UCHAR_MAX+1;    /* maximum value for a pixel, 255, plus 1 for rounding */
  128.     if(tmax<1.0)tmax=1.0;
  129.     ifunScale = LONG_MAX/tmax;
  130.     for (i=0;i<dim;i++) ifun[i] = f[i]*ifunScale;
  131.     
  132.     #if DIAGNOSTICS
  133.         /* This is a diagnostic printout to track down overflow problems */
  134.         t = 0L;
  135.         a = 0.0;
  136.         for (i=0;i<dim;i++) {
  137.             t += labs(UCHAR_MAX*ifun[i]);
  138.             a += fabs(f[i]);
  139.         }
  140.         device=GetGDevice();
  141.         SetGDevice(GetMainDevice());
  142.         PmForeColor(UCHAR_MAX);
  143.         PmBackColor(0);
  144.         PrintfExit("a=%lf, t=%lx, ifunScale=%lx, t/ifunScale=%lf\n",a,t,ifunScale,t*1.0/ifunScale);
  145.         SetGDevice(device);
  146.     #endif
  147.     
  148.     /*
  149.     Unfortunately CopyBits will crash if the line buffer is too large, e.g.
  150.     640 pixels long. We can minimize this problem by reducing the size of the
  151.     srcRect to be no larger than necessary (i.e. dstRect plus dim).
  152.     */
  153.     mysrcRect = *srcRectPtr;
  154.     mydstRect = *dstRectPtr;
  155.     if(convolveX){
  156.         i=mysrcRect.right-mysrcRect.left - (dim+mydstRect.right-mydstRect.left);
  157.         if(i>0)InsetRect(&mysrcRect,i/2,0);
  158.         i=mysrcRect.bottom-mysrcRect.top - (mydstRect.bottom-mydstRect.top);
  159.         if(i>0)InsetRect(&mysrcRect,0,i/2);
  160.         if(i<0)InsetRect(&mydstRect,0,-i/2);
  161.     }else{
  162.         i=mysrcRect.right-mysrcRect.left - (mydstRect.right-mydstRect.left);
  163.         if(i>0)InsetRect(&mysrcRect,i/2,0);
  164.         if(i<0)InsetRect(&mydstRect,-i/2,0);
  165.         i=mysrcRect.bottom-mysrcRect.top - (dim+mydstRect.bottom-mydstRect.top);
  166.         if(i>0)InsetRect(&mysrcRect,0,i/2);
  167.     }
  168.     srcRectPtr = &mysrcRect;
  169.     dstRectPtr = &mydstRect;
  170.  
  171.     if(convolveX){
  172.         dstLines = dstRectPtr->bottom-dstRectPtr->top;
  173.         srcLines = srcRectPtr->bottom-srcRectPtr->top;
  174.         dstLength = dstRectPtr->right-dstRectPtr->left;
  175.         srcLength = srcRectPtr->right-srcRectPtr->left;
  176.     }else{
  177.         dstLength = dstRectPtr->bottom-dstRectPtr->top;
  178.         srcLength = srcRectPtr->bottom-srcRectPtr->top;
  179.         dstLines = dstRectPtr->right-dstRectPtr->left;
  180.         srcLines = srcRectPtr->right-srcRectPtr->left;
  181.     }
  182.     if(dstLines > srcLines) lines = srcLines; 
  183.     else lines = dstLines;
  184.     
  185.     /* Allocate our line buffers. */
  186.     srcTmp = **(*GetGDevice())->gdPMap;    /* Use current device to init PixMap fields */
  187.     srcTmp.pmVersion=srcTmp.packType=srcTmp.packSize=0;
  188.     srcTmp.planeBytes=srcTmp.pmReserved=0;
  189.     srcTmp.cmpCount=1;
  190.     srcTmp.pixelSize=8;                    /* The main loop assumes 8 bits per pixel */
  191.     if(convolveX){
  192.         SetRect(&srcTmp.bounds,0,0,srcLength,1);    /* a horizontal line buffer */
  193.     }else{
  194.         SetRect(&srcTmp.bounds,0,0,1,srcLength);    /* a vertical line buffer */
  195.     }
  196.     srcTmp.rowBytes = 2*(((srcTmp.bounds.right-srcTmp.bounds.left)*srcTmp.pixelSize+15)/16);
  197.     srcTmp.baseAddr = (Ptr) malloc(srcTmp.rowBytes*(srcTmp.bounds.bottom-srcTmp.bounds.top));
  198.     srcTmp.rowBytes |= 0x8000;    /* Mark it as a PixMap, not a BitMap */
  199.     dstTmp = srcTmp;
  200.     if(convolveX){
  201.         SetRect(&dstTmp.bounds,0,0,dstLength,1);    /* a horizontal line buffer */
  202.     }else{
  203.         SetRect(&dstTmp.bounds,0,0,1,dstLength);    /* a vertical line buffer */
  204.     }
  205.     dstTmp.rowBytes = 2*(((dstTmp.bounds.right-dstTmp.bounds.left)*dstTmp.pixelSize+15)/16);
  206.     dstTmp.baseAddr = (Ptr) malloc(dstTmp.rowBytes*(dstTmp.bounds.bottom-dstTmp.bounds.top));
  207.     dstTmp.rowBytes |= 0x8000;    /* Mark it as a PixMap, not a BitMap */
  208.     if(srcTmp.baseAddr == NULL || dstTmp.baseAddr == NULL)
  209.         PrintfExit("ConvolveX/Y: Sorry, couldn't allocate enough memory for working arrays\007\n");
  210.     
  211.     /* Define the rectangles that map src and dst to their respective line buffers */
  212.     srcRect = srcTmp.bounds;
  213.     OffsetRect(&srcRect,srcRectPtr->left,srcRectPtr->top);
  214.     dstRect = dstTmp.bounds;
  215.     OffsetRect(&dstRect,dstRectPtr->left,dstRectPtr->top);
  216.     if(convolveX){
  217.         OffsetRect(&srcRect,0,srcLines/2-lines/2);
  218.         OffsetRect(&dstRect,0,dstLines/2-lines/2);
  219.     }else{
  220.         OffsetRect(&srcRect,srcLines/2-lines/2,0);
  221.         OffsetRect(&dstRect,dstLines/2-lines/2,0);
  222.     }
  223.     
  224.     /* Convolve one line per iteration */
  225.     /* y indicates which line, x indicates which pixel along that line */
  226.     for (y= -lines/2; y<(lines+1)/2; y++) {
  227.         /* first copy a line from srcBits into line-buffer srcTmp */
  228.         if(srcPixelSize == srcTmp.pixelSize)
  229.             CopyBitsQuickly(srcBits,(BitMap *) &srcTmp,&srcRect,&srcTmp.bounds,srcCopy,NULL);
  230.         else
  231.             CopyBits(srcBits,(BitMap *) &srcTmp,&srcRect,&srcTmp.bounds,srcCopy,NULL);
  232.         if(can32)SwapMMUMode((void *)&mmuMode);
  233.         if(!convolveX){
  234.             /* Unfortunately the minimum PixMap rowBytes is 2 bytes, so we
  235.             squeeze out the extra space before processing
  236.             */
  237.             srcPtr= (unsigned char *) srcTmp.baseAddr;
  238.             dstPtr= (unsigned char *) srcTmp.baseAddr;
  239.             for (x=0; x<srcLength; x++) {
  240.                 *dstPtr++ = *srcPtr;
  241.                 srcPtr += 2;
  242.             }
  243.         }
  244.         dstPtr= (unsigned char *) dstTmp.baseAddr;
  245.         srcx = srcLength/2 - dim/2 - dstLength/2;
  246.         for (x=0; x<dstLength; x++) {
  247.             imax=dim;
  248.             if(srcx+dim>srcLength) imax=srcLength-srcx;
  249.             if(srcx<0){
  250.                 imin= -srcx;
  251.                 srcPtr=(unsigned char *) srcTmp.baseAddr;
  252.             }
  253.             else {
  254.                 imin=0;
  255.                 srcPtr=(unsigned char *) srcTmp.baseAddr+srcx;
  256.             }
  257.             fPtr=&ifun[imin];
  258.             for (t=0L, i=imin;i<imax;i++) t += *fPtr++ * (long)(unsigned long) *srcPtr++;
  259. #if IMPROVED_ROUNDING
  260.             if(t>0)t+=ifunScale>>1;        // round t/ifunScale to nearest integer
  261.             else t-=ifunScale>>1;
  262. #endif
  263.             t/=ifunScale;
  264.             if(convolveX){
  265.                 *dstPtr++ = (unsigned char) t;
  266.             }else{
  267.                 /* Unfortunately the minimum PixMap rowBytes is 2 bytes, so we
  268.                 insert a space between pixels
  269.                 */
  270.                 *dstPtr = (unsigned char) t;
  271.                 dstPtr += 2;
  272.             }
  273.             srcx++;
  274.         }
  275.         if(can32)SwapMMUMode((void *)&mmuMode);
  276.         /* finally, copy the line from buffer dstTmp to dstBits */
  277.         if(dstPixelSize == dstTmp.pixelSize)
  278.             CopyBitsQuickly((BitMap *) &dstTmp,dstBits,&dstTmp.bounds,&dstRect,srcCopy,NULL);
  279.         else
  280.             CopyBits((BitMap *) &dstTmp,dstBits,&dstTmp.bounds,&dstRect,srcCopy,NULL);
  281.         if(convolveX){
  282.             OffsetRect(&srcRect,0,1);    /* adjust Rects to point to next line */
  283.             OffsetRect(&dstRect,0,1);
  284.         }else{
  285.             OffsetRect(&srcRect,1,0);    /* adjust Rects to point to next line */
  286.             OffsetRect(&dstRect,1,0);
  287.         }
  288.     }
  289.     free((void *) ifun);
  290.     free((void *) srcTmp.baseAddr);
  291.     free((void *) dstTmp.baseAddr);
  292. }
  293.  
  294.